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Abstract 

In this numerical study, we investigate the role of intrinsic heterogeneities of cardiac tissue due to M cells 
in the generation and maintenance of reentrant excitations using the detailed Luo-Rudy dynamic model. 
This model has been extended to include a description of the long QT 3 syndrome, and is studied in both 
one dimension, corresponding to a cable traversing the ventricular wall, and two dimensions, representing 
a transmural slice. We focus on two possible mechanisms for the generation of reentrant events. We first 
investigate if early-after-depolarizations occurring in M cells can initiate reentry. We find that, even for 
large values of the long QT strength, the electrotonic coupling between neighboring cells prevents early- 
after-depolarizations from creating a reentry. We then study whether M cell domains, with their slow 
repolarization, can function as wave blocks for premature stimuli. We find that the inclusion of an M cell 
domain can result in some cases in reentrant excitations and we determine the lifetime of the reentry as a 
function of the size and geometry of the domain and of the strength of the long QT syndrome. 
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Spatial heterogeneity of cardiac tissue causes cells to repolarize at different rates, leading to 
a dispersion of repolarization time. Dispersion of repolarization has been postulated to function 
as a substrate for reentry phenomena which occur when the propagation of the electric wave 
is blocked in one direction causing the wave front to curl and reenter the previously excited 
tissue. The ensuing incoherent electrical activity of the heart is thought to subsequently lead 
to ventricular fibrillation, the leading cause of sudden death in the industrialized world. In 
this numerical study, we investigate the role of intrinsic heterogeneities due to M cells in the 
generation and maintenance of reentrant excitations. These cells, located in the midmyocardium, 
exhibit a prolonged action potential duration and the resulting dispersion of repolarization is 
exacerbated in patients with the long QT syndrome. Using a detailed electrophysiological model, 
two previously postulated mechanisms are investigated: one in which the M cells create premature 
stimuli that lead to reentrant events and one in which the M cell domains function as wave blocks 
for electrically propagating waves originating, as usually, from the endocardium. We find that 
the first mechanism is unlikely to happen under the conditions we study, due to the strong 
electrotonic coupling that is normally present between neighboring cells. The second mechanism, 
on the other hand, can result in reentrant excitations. We determine the lifetime of the reentry as 
a function of the size and geometry of the domain, and of the strength of the long QT syndrome. 

1 Introduction 

Spatial heterogeneity of cardiac cells is commonly believed to play a major role in the generation and maintenance 
of arrhythmia. Both dynamic heterogeneity, where the dynamical properties of the homogeneous tissue are 
responsible for the heterogeneity, 1-3 and intrinsic heterogeneity, due, for instance, to the presence of different 
types of cells in the cardiac tissue, can result in cardiac cells that repolarize at different rates. The resulting 
dispersion of repolarization can lead to regions that block propagating waves and can, ultimately, result in 
reentry phenomena. 

The focus of this study is intrinsic heterogeneities, which are present in a variety of forms. One of the most 
dramatic intrinsic heterogeneities is the presence of a layer of cells, M cells, roughly located in the middle of the 
myocardial wall. 4 ' 5 This layer consists of cells that have distinct electrophysiological properties, most notably a 
prolonged action potential duration. The presence of these cells has been observed in many experimental studies 
and has been found in the ventricles of the dog, 4 guinea pig 6 and human heart. 7 M cells may constitute up to 
40% of the human left ventricular free wall 7 and the resulting action potential prolongation is most dramatic 
for isolated cells but remains substantial in wedge preparations where the M-cells are coupled through gap 
junctions to the cells in the epi- and endocardium. 8 The dispersion becomes even more dramatic under the 
influence of cardiac drugs that prolong the action potential and/or in the presence of slow pacing. 5 The ion 
channels responsible for the prolongation of the action potential in M cells include primarily a smaller slow 
activating potassium current Ik s , an augmented late sodium current I Na and a larger sodium-calcium exchange 
current. 9 

It has been hypothesized that the resulting transmural dispersion plays a major role in the initiation and 
maintenance of reentry. 10 In particular, the role of M cells in the initiation of torsade de pointes, a polymorphic 
ventricular tachycardia, in patients with the long QT syndrome has recently received considerable attention. 11 
The long QT syndrome is characterized by a prolonged ventricular repolarization and is associated with a high 
risk of sudden cardiac death. It manifests itself as an abnormally long QT interval in the ECG and can be either 
congenital, idiopathic or acquired (for a review see Ref. 12). Several forms of long QT syndromes with different 
cellular mechanisms and clinical characteristics have been identified. Nearly all cardiac events in patients with 
long QT 1, associated with a reduced Ik s , are initiated during increased sympathetic activation induced by 
physical or emotional stress. 13 ' 14 On the other hand, in patients with long QT 2, characterized by a decrease 
in Ixr, and long QT 3, corresponding to a late In<i, cardiac events occur most often during rest or sleep. 1 

Several possible mechanisms linking M cells and dispersion of repolarization with cardiac arrhythmias have 
been postulated. In one, M cells develop early-after-depolarizations, which are characterized by a depolarization 
at the end of the action potential plateau. This depolarization could then create a premature stimulus via the 
electrotonic interaction with neighboring cells. This premature stimulus will be partially blocked by the M cell 
domain but can initiate a wave in other directions. Once the repolarization is complete, this wave can reenter 
the M cell zone, leading to a classic figure-of-eight reentry. 15 Crucial to this mechanism is the ability for the 
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early-after-depolarizations to excite neighboring cells despite electrotonic effects that will tend to smooth out 
steep gradients. 

In another possible mechanism, the dispersion of repolarization due to M cells creates regions with delayed 
repolarization that act as conduction blocks for subsequent waves. Under suitable conditions, the dispersion 
of repolarization can then lead to partial wave block and reentry. This scenario has recently been shown to 
be plausible in a canine wedge preparation of the left ventricle. 11 In this study, a transmural cross section of 
the wedge was visualized via voltage sensitive dyes and long QT 2 conditions were realized pharmacologically 
by adding the Jr> blocker d-sotalol. 16 This led to a marked increase in dispersion of repolarization and the 
presence of distinct islands of M cells, which resulted in areas with steep spatial gradients of repolarization. 
When the wedge was paced under bradycardic conditions from the endocardium, followed by a premature 
stimulus delivered at a short coupling interval (time interval between the premature stimulus and the previous, 
regular stimulus), the islands were found to function as conduction blocks and torsade de pointes ensued. 
The underlying mechanism of the resulting arrhythmia was demonstrated to be sustained intramural reentrant 
excitations. 

Only a limited number of numerical studies investigating long QT syndromes, dispersion of repolarizations 
and its role in arrhythmias has been undertaken. Work by Rudy and colleagues investigated long QT syndrome 
in isolated cells 17 ' 18 and strands of cells with embedded domains of M cells 19,20 but did not address the role 
of M cells in the initiation of reentry. Clayton et al. 21 studied the size of the vulnerable window, defined 
as the time interval for which a premature stimulus initiates a unidirectional wave, in homogeneous strands 
of tissue. They found that long QT conditions did not increase this size measurably. However, they did not 
investigate propagation block in strands having domains consisting of electrophysiologically different cells. They 
also investigated the tip trajectory in numerical models for long-QT 1, 2 and 3, using a relatively large square 
homogeneous domain but did not study the initiation of reentrant excitation. 

In this paper, we use a detailed electrophysiological model 22 to investigate the two possible mechanisms 
detailed above. We first study if normally coupled tissue can generate: unidirectional waves in one-dimensional 
(ID) strands of tissue. We then determine under what conditions a domain of M cells will block propagation 
in a ID cable. We finally determine under what conditions a two-dimensional (2D) sheet of tissue containing a 
layer of M cells will lead to reentrant excitation. 



2 The Model 

To describe the electrophysiological properties of the cardiac cell we use the Luo-Rudy dynamic (LRd) model. 22, 23 
The dynamics of the tissue is described by the reaction-diffusion equation 

— = — - — W 2 V-— (1) 

Ot C m pS v Cm 

where V (mV) is the membrane potential, C m (/zFcm~ 2 ) is the membrane capacitance, p (ficm) is the bulk 
resistivity, S v (cm -1 ) is the surface to volume ratio and V 2 is the Laplacian operator. The combination of 
c ^ g gives the diffusion constant D, which we have taken to be D — c x pS =0.11 cm 2 /s resulting in a planar 
wave speed of approximately 17 cm/s. This speed is consistent with experimentally observed wavespeeds in 
canine ventricles when measured perpendicular to fibers. 24 As is the case for all electrophysiological models, 
the details of the LRd model are contained in the description of the total ionic current Ii on (/iAcm -2 ) flowing 
through the cardiac cell membrane. In contrast to earlier models, which describe only a small number of ionic 
currents, the LRd model and other often so-called "second generation" models attempt to incorporate as many 
significant currents as possible, including a more complete description of intracellular calcium dynamics. 

Our choice of the LRd model was motivated by several considerations. First, all other existing models we 
tested, including the first generation Luo-Rudy 25 and Beelcr-Rcutcr 26 models and the second generation model 
of Fox et al., 27 exhibited spiral cores that were typically comparable to or larger than the computational domain, 
taken to be typical of human ventricles. Thus, these models are not suitable to investigate sustained reentry in 
our computations. Second, we are interested in modifying specific currents that are responsible for the distinct 
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Figure 1: (a) Typical profile of the maximal conductance of Iks, Gks, used in the numerical studies of a cable. 
The conductance is shown normalized by the maximal conductance in normal tissue, GksO- The cable consists 
of a region of width w with a reduced Gks (H) corresponding to M cells, sandwiched between two equally 
sized domains of normal tissue (I and III), (b) The geometry considered in our 2D simulations, consisting of a 
region of M cells, shaded in gray, surrounded by normal tissue. As in the cable, the M cells are modeled via the 
reduction of Gk s and the transition between the different domains is taken to be smooth. 



properties of M cells and that play a role in the long QT syndrome (see below). The LRd model, being a second 
generation model, allows us to change these currents directly. 

Since M cells are predominantly present in the midmyocardium, we focus here on the effects of regional 
heterogeneities in transmural geometries. Furthermore, since the transmural wedge experiment in Ref. 11 
suggests that, at least in the initial stages, the reentry caused by the long QT syndrome is a 2D phenomenon 
we have limited ourselves to tissue sheets that represent thin transmural wedges. We have taken the transmural 
dimension to be 1 cm, while the other dimension is taken to be 2.5 cm. Also, we found it useful to investigate 
the behavior of the tissue in a ID cable that can be thought of as a straight line traversing the myocardial wall 
and perpendicular to the endo and epicardial walls. In all simulations, we have taken a spatial discretization 
of 0.008 cm and a time step of 0.02 ms. We have verified that using smaller time and space steps did not 
significantly modify the quantitative results. 

Spatial heterogeneity was introduced by subdividing the cable into three domains (see Fig. QJi). Since we 
focus here primarily on the role of the M cell domain, we have chosen two outer domains that have equal length 
and identical electrophysiological properties. This is, of course, a simplification as the endo- and epicardium can 
have different electrophysiological properties. 9 However, in the canine wedge preparation, the action potential 
durations of the epi- and endocardial cells were only minimally different. 11 The middle domain consists of a 
region of M cells, which, electrophysiological^, were modeled by decreasing the value of the maximum conduc- 
tance Gks of the slow activating potassium current Iks- We tuned the value of Gks such that the epi- and 
endocardial cells and the M cells had an action potential that correspond to the experimental findings in Ref. 
11. Specifically, we took the normal tissue to have a maximal conductance GksO that was a factor of 1.5 larger 
than the original value in the LRd model while the maximal conductance GksI in the M cells was chosen as 
GksI — 0.375GksO- 

The transition between the two regions was made smooth by employing a simple sigmoidal function for 
Gks (see Fig. QJi and supplemental material). The two-dimensional (2D) geometry consists of a transmural 
slice of the myocardial wall containing a strip of M cells. Of course, it is impossible to initiate reentry via 
uniform pacing of the entire endocardium if the computational domain is invariant under translation along the 
endocardium (the vertical direction in Fig. \2p). Thus, to break the symmetry we have chosen to use an M cell 
layer consisting of a strip of width w that starts at one boundary and ends with a "bubble" with radius r at the 
other end (see Fig. QJ)). Evidence for inhomogeneous distributions of M cells can be seen in ref. 11 where the 
repolarization maps show distinct islands of slower repolarization, presumably corresponding to M cell domains. 
As in the ID cable, the transition between the different domains was made smooth. 

We will focus here mostly on long QT 3, which is characterized by an incomplete inactivation of the sodium 
current I^a and which leads to an abnormally prolonged action potential, particularly in M cells. To model the 
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Ijsr a channel and to allow a description of the long QT 3 mutant we have adopted the strategy of Clancy and 
Rudy and have replaced I^a in the LRd model with a Markovian description. 28 That is, instead of calculating 
the opening and closing of the fast sodium channel using gating variables that are deterministic functions of 
the membrane potential, we describe this channel using a collection of open, closed and inactivation states. 
The transition rates between these states are dynamically calculated and are functions of V. In addition to 
the background mode, also present in wild type cells, mutant cells have a burst mode which does not include 
inactivation states. When the channel is in this burst mode, the inactivation transiently fails, leading to a leak 
current and a prolonged action potential. 

To simplify the original model for Ino 28 and its subsequent refinement, 29 and to render our computations 
more efficient, we have changed the original Markovian model. More details of this modification can be found in 
the supplemental material. In essence, the simplification uses the existence of two distinct timescales of the I^ a 
channel: a short timescale on the order of a single action potential and a much longer timescale on the order 
of multiple action potentials. Using these two timescales the Markov model can be reduced to an equation for 
Inu which is similar to the original formulation in the LRd model supplemented by a dynamical equation which 
describes the fraction of channels in the burst mode. Numerical checks, presented in the supplemental material, 
show that the reduced description accurately reproduces the leak current. Furthermore, we have verified that 
action potentials generated with the simplified model do not differ significantly from those generated using the 
full model. 

The simplification of Ino,, along with the geometry of our computational domain, lead to several adjustable 
parameters in our simulations. First, the strength of the long QT syndrome can be changed via a single 
dimensionless parameter /i. It is a measure of the fraction of mutant cells (see appendix) and is equal to zero 
in the case of normal cells. Second, the width of the M cell layer w can be changed. Finally, in 2D we have the 
additional parameter r, the radius of the disk at the end of the M cell layer. Notice that the minimal value of 
r, r = w/2, corresponds to a strip with a rounded edge. 

3 Results 

One dimension 

We first quantified the effect of the size of the M cell domain w and of the long QT syndrome strength \i on 
action potential duration in a cable. For this, we stimulated the cable at one end with a fixed basic cycle length 
of 1000 ms for 40 s. We then measured the action potential duration in the cell at the midpoint of the cable 
(and thus at the midpoint of the M cell domain) and at the end of the cable. In Fig. [2^ we plot the action 
potential duration as a function of \x at the midpoint of a cable with (solid line, w = 0.48 cm) and without 
an M cell domain (dashed line, w = 0.0 cm). Fig. |2h shows the corresponding repolarization times along the 
cable, measured as the time interval between the stimulus and the time at which the membrane potential has 
repolarized to V — —70 mV. As expected, increasing the long QT syndrome strength results in an increase in 
action potential durations. Since the action potential duration in M cells is more strongly affected by an increase 
in //, higher values of fi result in an increased dispersion of repolarization time in the tissue. In Fig. we 
show the action potential durations at the midpoint of the heterogeneous cable as a function of w for fi — 0.08, 
the long QT syndrome strength employed in most of this study. The corresponding repolarization times can 
be found in Fig. |2Jl. The action potential duration of points well into the normal tissue are hardly affected by 
increasing fi while the midpoint of the cable displays a marked increase in action potential duration when w is 
increased. The action potential durations computed numerically here for /j, = and fj, = 0.08 are similar to the 
ones observed in Ref. 11 when the wedge was paced under bradycardic conditions with and without d-sotalol, 
respectively. 

The shape of the action potential at different sites of the cable can be seen in Fig. [3] Here, we have divided 
the cable into 7 subdomains and have plotted the action potential for each subdomain for a fixed stimulation 
period of 1000 ms and varying long QT syndrome strength. The domain of M cells has a width of w — 0.48 
cm and corresponds to the region c, d and e. It can clearly be seen that the long QT syndrome has a more 
dramatic effect on the M cells (d) than on normal cells (a) and (g) . Also note that due to the symmetry breaking 
introduced by the pacing at the left end of the cable, the action potentials are longer at the right boundary of 
the M cell domain (f) than at the left boundary (b). 

Next, we investigated if the presence of M cells could provoke early after depolarizations and waves prop- 
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Figure 2: (a): Action potential duration at a basic cycle length of 1000 ms as a function of the strength of the 
leak //va, (i, at the midpoint of a cable with (solid line, w = 0.48 cm) and without an M cell domain (dashed 
line), (b): Corresponding repolarization times for different values of /j, varying from (bottom line) to 0.09 
with steps of 0.01. (c): Action potential duration at the midpoint of the cable as a function of the width of the 
M cell domain (basic cycle length=1000 ms and /i = 0.08). (d): Corresponding repolarization time for different 
widths of the M cells region ranging from a homogeneous cable (w = cm and bottom line) to a domain of 
width w = 0.48 cm with increments of 0.06 cm. 




Figure 3: Action potentials at different positions along a cable of length 1 cm for increasing values of /j, from 
(shortest action potential) to 0.09 (longest action potential) with steps of 0.01. The domain of M cells has a 
width of w = 0.48 cm and corresponds to the region c, d and e, while the transition regions are b and f. The 
basic cycle length is 1000 ms, the horizontal bar corresponds to 200 ms and the vertical bar corresponds to 50 
mV. 



6 



agating unidirectionally in our cable. We first verified that single M cells exhibit early after depolarizations, 
defined as an increase in the membrane potential during the plateau phase of the action potential. We found 
that these early after depolarizations are present for a relatively wide range of parameters and pacing protocols 
and have an amplitude that is consistent with experimental findings. This is shown in Fig. 0^,1 where we plot, as 
a solid line, three subsequent action potentials for an isolated M cell with /i =0.06. In contrast, an isolated epi- 
or endocardial cell, paced at the same basic cycle length of 1000 ms, does not exhibit early after depolarizations 
(dashed line). 

Having established the presence of early after depolarizations in isolated cells, we next considered a cable 
containing a domain of M cells. In Fig. 0Jl2, we show the action potential in a cable consisting of cells identical 
to the ones of Fig. The solid line represents the membrane potential at the middle of the M cell domain 

(w — 0.48 cm), while the dashed line represents the membrane potential in the middle of the epicardial domain. 
The early after depolarizations, present in the isolated cells, have now disappeared due to electrotonic effects. 
We found qualitatively similar results for different basic cycle lengths, long QT strengths and widths of the 
M cell domain: early after depolarizations provoked in single M cells were always found to be suppressed by 
electrotonic effects once the cells were inserted in a cable containing an M cell domain sandwiched between an 
endo- and epicardial domain. 

In fact, we found that, for the model used here, the only possible way to generate an early after depolarization 
that travels along a homogeneously coupled cable is to either change the pacing protocol or to increase the long 
QT strength significantly. Both scenarios lead to early after depolarizations that are already present in the 
endocardial domain. An example is shown in In Fig. 0Jd where we have plotted the action potential at four 
different locations along the cable as indicated in the figure. The cable was paced at a constant basic cycle length 
of 1000 ms, followed by a stimulus with coupling interval of 2000 ms. The endocardial domain now generates 
an early after depolarization with a significant amplitude that is able to travel through the M cell domain into 
the epicardial domain. However, the early after depolarization is propagating with a smaller speed that the 
repolarization wave and disappears before reaching the end of the cable. Thus, in conclusion, under normal 
coupling conditions, the diffusive term in Eq. JTJ always leads to considerable dissipation of the membrane 
potential, which prevents the early after depolarizations from stimulating neighboring tissue. 

Some studies have suggested that a discontinuity in conductivity exists between the M cell domain and the 
epicardium, leading to inhomogeneous coupling. 8 It is possible that the inclusion of such inhomogeneous cou- 
pling, as shown in other modeling studies that investigated the propagation from the Purkinje fibers (specialized 
fibers that rapidly transmit impulses from the atrioventricular node to the ventricles) to the mycardium, 30 ' 31 al- 
lows for unidirectional conduction. However, we found that propagation provoked by early after depolarizations 
in our model was only possible when we altered the conductance significantly. Specifically, we found that we had 
to reduce the conductance by at least a factor of ten, which is much more that the factor of three reported in 
Ref. 8. An example is shown in Fig. [SJwhere we have plotted the membrane potential at four different locations 
in the cable, as indicated in the figure. An early after depolarization appears in the M cell domain (b) as the 
coupling interval is increased. This early after depolarization becomes more pronounced at the boundary of 
the M cell domain and the epicardium, where the resistance is increased by a factor of thirteen. The increased 
resistance allows the epicardium to repolarize and the early after depolarization eventually creates an extra 
stimulus in the epicardium. 

These results suggest that the first arrhythmic mechanism described in the Introduction is unlikely to occur 
in the model we consider here. To study the second mechanism, we recall that this scenario predicts the blocking 
of a wavefront by the remnants of previous stimuli. To investigate the possibility of a wave block created by 
the M cell domain, we again paced the cable at one end with a fixed period. After 40 cycles, we delivered a 
premature stimulus to the same end of the cable and varied the coupling interval between the last of the regular 
cycles and the premature stimulus. Conduction block was defined as an inability to propagate through the M 
cell domain to the opposite end of the cable. The experiment was repeated for different basic cycle lengths 
and the conduction block window, the range of coupling intervals that resulted in a conduction block, was 
determined. 

We found a rich set of possible dynamics that is summarized in Fig. where we have plotted our findings 
in the basic cycle length (BCL) vs. strength of the long QT syndrome (fi) phase space. For wild type cells, 
i.e. fj, — 0, the presence of the M cells domain does not lead to any conduction block. However, qualitatively 
different wave propagations are present. For small pacing periods, the premature wave travels through the M 
cell region with almost unchanged speed. This behavior is demonstrated in Fig. EJl which shows a color-coded 
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Figure 4: al: The membrane potential for a single isolated M cell (solid line) for a long QT strength /j, — 0.06, 
along with the membrane potential for a single isolated epi- or endocardial cell (dashed line). The cell is 
paced with a basic cycle length equal to 1000 ms. a2: Membrane potential at different positions along our 
inhomogeneous cable (w=0A8 cm) consisting of cells used in al. The solid line correspond to the midpoint 
of the M cell domain, while the dashed line correspond to the middle of the epicardial domain. bl-b4: The 
membrane potential at four different locations along the cable, as indicated by the graph on the right. The M 
cell domain (w=0A8 cm) is shaded gray and the cable (/z=0.09) is paced using a basic cycle length of 1000 ms, 
followed by a stimulus with a coupling interval of 2000 ms. 




Figure 5: Action potentials at different positions along a cable of length 1 cm shown on the right. The shaded 
M cell domain has width w — 0.6 cm and the cable (p = 0.08) is paced with basic cycle length of 1000 ms, 
followed by a stimulus with coupling interval of 1200 ms. The maximal conductance Gk s of the slow activating 
potassium current was increased by a factor of two in the epicardial domain. The resistance between the last 
gridpoint of the M cell domain and the first gridpoint of the epicardium was abruptly increased by a factor of 
thirteen and returned to normal over the following ten gridpoints. 
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Figure 6: Phase diagram of the response of a cable paced at a regular basic cycle length (BCL) to a premature 
stimulus in the presence of a M cell domain of width 0.48 cm and of varying long QT syndrome strength. In 
regions I and I* no conduction block occurs: a wave that has begun to propagate will propagate through the 
cable. In region I the wave speed in the M cell domain is similar to the wave speed in the normal domain while 
in region I* the wave in the M cell domain slows down considerably before resuming normal propagation. The 
dashed line is a guide to the eye to separate the two different regions. Regions II and III correspond to two 
different types of conduction block. In region II the normally propagating wave reaches the M cell region and 
fails to propagate almost immediately. Region III, on the other hand, corresponds to conduction block in which 
the wave is able to enter the M cell domain, where it propagates with a significantly decreased wave speed before 
failing to propagate. The size of the window of coupling intervals which lead to conduction block in region II 
and III is indicated by a gray scale, with white being the smallest range and black being the largest range. 



space-time plot of the membrane voltage. The parameter values for Fig. |7Ji correspond to the point marked 
a in Fig. H3 Its behavior is typical of all points in region I in Fig. H3 For longer pacing periods, the wave 
front slows down significantly when entering the M cell region. After coming almost completely to a stop, with 
typical wave speeds found to be around 4 cm/s, the wave resumes its original speed and propagates into the 
epicardium. This slow- down- and- go behavior is demonstrated in Fig. 0:, corresponding to point c in Fig. H3 
and is representative of the points in the region labeled I* in Fig. [Gj 

For sufficiently large values of fi, conduction block occurs (region II and III in Fig. EJ). Two qualitatively 
different types of conduction block are observed. In region II, the wave penetrates the M cell domain only 
slightly before being blocked. A space-time plot of point b within region II is shown in Fig. 0d. The conduction 
block window for points in region II is rather small, typically less than 20 ms. This range is indicated in Fig. 
on a linear gray scale with white corresponding to small conduction block windows and black corresponding to 
large ones (> 30 ms). The conduction block window for point b was found to be 15 ms. 

The second type of conduction block can be found in region III and is shown in Fig. 0i. A wave front 
enters the M cell domain, reduces its speed significantly and is subsequently blocked. The conduction block 
window for this region is much larger than for region II. For example, the conduction block window for point d 
is 72 ms. This is also shown in Fig. [S] where we have plotted the conduction block window as a function of the 
width of the M cell domain for two different values of the long QT strength and for two different values of the 
pacing period. For both regions II and III the conduction block windows increases as w is increased. We have 
also investigated the location of the boundaries of the regions in Fig. for different w. Increasing w did not 
significantly increase the size of any of the regions, while decreasing w reduced the size of region II and III. Of 
course, for w = no conduction block occurs and both region II and III are absent in the phase diagram. 

In most of region III the behavior of the wave as a function of the coupling interval for fixed basic cycle 
length and \i can be summarized as follows: For large coupling intervals, the wave is able to propagate normally 
through the M cell domain. Upon decreasing the coupling interval the wave slows down significantly inside the 
M cell domain. Below a critical coupling interval conduction stops within the M cell domain. Finally, after 



9 



50 mV 
OmV 
-50 mV 

-100 mV 

x 

Figure 7: Space-time diagrams of the transmembrane potential of a cable of length 1 cm that includes an M cell 
domain of width w = 0.48 cm and long QT strength fi — 0.07. The cable is paced at the left end with constant 
basic cycle length for 10 s, after which a premature stimulus of varying coupling intervals is given. The color 
code used is shown by the colorbar on the right. The four figures a-d correspond to the points marked a-d in 
Fig. Eland are typical of the dynamics in the four different regions in that figure, (a) Dynamics representative 
of region I in Fig. [5] the premature stimulus (coupling interval=140 ms) propagates through the M cell domain 
with almost the same speed (basic cycle length=300 ms). (b) Dynamics Representative of region II in Fig. 
the wave is blocked immediately when entered the M cell domain (basic cycle length=500 ms and coupling 
interval=184 ms). For this basic cycle length, the conduction block window was fairly small (15 ms) with 
conduction block occurring for a coupling interval between 170 ms and 185 ms. (c) Dynamics Representative 
of region I* in Fig. H3 due to the slow repolarization of the M cells the wave speed is significantly decreased 
when entering the M cell domain. The speed increases again when the wave can propagate freely through the 
repolarized medium (basic cycle length=800 ms and coupling interval=224 ms). (d) Dynamics Representative 
of region III in Fig. |SJ the wave speed decreases dramatically upon entering the M cell domain. However, unlike 
(c), normal propagation is unable to resume and the wave is blocked (basic cycle length=1000 ms and coupling 
interval= 245 ms). For this stimulation period, we found a conduction block window of 72 ms with conduction 
block occurring for a coupling interval between 225 ms and 297 ms. 
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Figure 8: The conduction block window duration (CBW) as a function of the width of the M cell domain for 
/i = 0.07 (dashed lines) and \x — 0.08 (solid lines) for a basic cycle length of 1000 ms (a) and 500 ms (b). Note 
the different scales of the CBW axis. 



10 



/ \ (cm) 

032- ■■■!■■] 

□ ■□□□J 

0.24- □□□□] 

□ □□J 
0.16- □□ 

0.08- n 

^ — I 1 h 




I 



■i— > 



2s 
1.5 s 
1 s 
0.5 s 
0s 



0.16 0.32 0.48 0.64 w ( cm ) 



Figure 9: Phase diagram: Duration of the reentry as a function of the two parameters governing the geometry 
of the 2D sheet: the width of the M cell region (w) and the radius of the bubble (r). The gray-scale indicates the 
duration as shown by the bar on the right (units are in seconds). A relatively sharp transition occurs between 
a regime for which reentries last at most around 1 second and a regime for which reentries last for more than 
two seconds and are pinned by the bubble (black squares). 



decreasing the coupling interval even further, propagation fails already in the epicardial domain (domain I in 
Fig. [IJi). Close to the boundary of region III with region I* a different behavior is observed. Again, for large 
coupling intervals the wave is able to propagate normally through the M cell domain. For smaller coupling 
intervals the wave slows down, and below a first critical coupling interval, propagation fails. Upon further 
decreasing the coupling interval, however, wave propagation resumes, characterized by a slowly propagating 
front within the M cell domain. Below a second critical coupling interval propagation fails in the epicardial 
domain. We will leave the discussion of this behavior to future studies. 

In summary, our ID results show that, for properly timed premature stimuli and sufficiently large values of 
H, the presence of a domain of M cells can create a conduction block. Furthermore, early-after-depolarizations 
were observed to develop within the M cell domain. However, in our numerical experiments, they are unable to 
excite already repolarized neighboring tissue due to electrotonic interactions. Thus, reentrant excitations will 
not be generated in our model via the first mechanism and in our 2D simulations we will exclusively focus on 
the second mechanism. 



Two dimensions 

To investigate whether the second mechanism could lead to reentrant excitation in 2D we performed an extensive 
numerical study of a transmural sheet, shown in Fig. ^p. Our simulation protocol was chosen to reduce the 
computational cost while minimizing any drift in the action potential duration. For this, we first divide the 2D 
sheet into ID cables running from the endo- to epicardium. These cables are then paced from the endocardial end 
for 40 s, after which the sheet is reassembled. Then, the entire sheet is paced for a further 10 s by stimulating the 
entire endocardium. This is finally followed by a premature stimulus with a varying coupling interval. Rather 
than choosing the premature stimulus at the epicardium as was done in Ref . 1 1 , we used a stimulus originating 
from the entire endocardium. This can be thought of as a premature stimulus coming from the Purkinje fiber. 31 
Our main result is summarized in Fig. where we show the length of reentry, measured in seconds, as a 
function of the two geometry parameters, w and r, for fixed long QT syndrome strength fi = 0.08. Of course, 
the duration of reentry is also dependent on the coupling interval (see below and Fig. I12f> and to determine its 
value, we varied the coupling interval in steps of 1 ms and measured the duration of reentry for the five coupling 
intervals that produced the longest lived reentrant excitations. The average of these five events is shown in 
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t(ms) 

Figure 10: Spiral tip trajectories for w = 0.48 cm and r = 0.24 cm (a) and w = 0.24 cm, r = 0.32 cm (b). The 
M cell domain is lightly shaded and the long QT syndrome strength is n = 0.08. In both cases, the start of 
the tip trajectory is indicated by an arrow and the initial stage of the reentrant excitation is not shown. Note 
that only a part of the computational domain is shown. In (a) , the dashed line is used as a guide to the eye 
in order to show that the spiral ends up leaving the sheet. The size of the tissue sheet is 1 cm x 2.5 cm. (c) 
Simulated pseudo-ECG corresponding to the reentry shown in (b), showing the 4th and 5th regular beat and 
the premature stimulus indicated by the arrow. 



Fig. [5] using a gray scale with white corresponding to short-lived reentry and black corresponding to long-lived 
reentry. Note that the diagonal corresponds to a strip with a rounded edge. 

From the figure, one can see that reentry that lasts longer than 1 s requires a sufficiently large M cell domain. 
This can be accomplished by either having a wide strip (w > 0.48 cm), or having a large bubble (r > 0.28 cm 
for w = 0.24 cm). These two cases, however, result in qualitatively different reentrant patterns. A typical spiral 
tip trajectory for a wide, rounded-edged, strip is shown in Fig. ITUk . The tip was calculated every 1 ms using 
the intersection of the iso-contour V = —30 mV and the line defined by h = 0.5, where h is the inactivation gate 
of the sodium channel. 21 The initial phase of the reentry is not shown and the first point of the tip trajectory 
is marked by the arrow. The spiral eventually drifts out of the computational domain, which is indicated by 
the dashed line. The total trajectory in Fig. HOfa represents 1 s and approximately 5 spiral rotations. 

The tip trajectory in the case of a large bubble is shown in Fig. 110b where we again have omitted the initial 
part of the reentrant excitation and marked the first point of the trajectory with an arrow. Here, the spiral 
tip becomes pinned to the bubble, leading to long lasting-reentry. A qualitative similar picture was found in 
a numerical study of a tissue sheet with a central ischemic region. 32 The spiral tip is plotted for roughly 2 s 
and 10 spiral rotations in Fig. 110b . The corresponding simulated pseudo-ECG, using the algorithm in, 33 is 
shown in Fig. 110b . In this figure, we have included the two regular preceding beats and have indicated the 
premature stimulus by an arrow. The polymorphic undulating ECG is suggestive of torsade de pointes. Finally, 
the reentrant excitation of Fig. 1 10b is also plotted in Fig. El where we show the activation of the tissue every 
25 ms using the color scale of Fig. 

Since reentry is caused by the partial blocking of the premature stimulus, we found a good correlation 
between the range of coupling intervals that leads to reentry and the range of coupling intervals that leads to a 
ID conduction block. This is illustrated in Fig. 1 121 where we have plotted the length of reentry as a function of 
the coupling interval for fi = 0.08, w = 0.32 cm and three different values of r. For these parameter values, the 
ID cable simulations reveal a conduction block for coupling intervals between 240 ms and 300 ms. In 2D, the 
onset of reentry occurs for a coupling interval of 240 ms, as can be seen in Fig. Eland no reentry is observed 
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Figure 11: Color plots of the membrane potential of the reentrant excitation of Fig. 110b . The M cell domain 
is indicated by the white dashed line. The plots are taken 25 ms apart and start Is after the application of the 
premature stimulus. The first one is in the upper left hand corner while the last one is in the lower right hand 
corner. The color scale is identical to the one employed in Fig. [7| Parameter values: /i=0.08, w = 0.24 cm, 
r = 0.32 cm and a physical domain of 1 cm x 2.5 cm 
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Figure 12: The length of reentry as a function of coupling interval for the transmural sheet with /j, — 0.08, basic 
cycle length=1000 ms, w = 0.32 cm and three different values of r (0.32cm (solid), 0.24cm (short dashed) and 
0.16cm (long dashed)). Reentrant excitations for r = 0.32 cm (solid line) that were characterized by a pinning 
to the bubble and that persisted throughout the duration of our numerical simulations (2 s) are indicated by 
the dotted line. 



for coupling intervals larger than 300ms. As already mentioned before, a larger bubble size leads to a more 
sustained reentry. In fact, for large bubbles, the spiral tip becomes attached to the bubble, resulting in reentry 
persisting throughout the full duration of our numerical simulations (2 s). These points are indicated by the 
dotted line in Fig. El 

4 Discussion 

We have presented a numerical study of the role of a particular type of spatial heterogeneity in cardiac tissue 
in the initiation and maintenance of reentrant excitations. Specifically, we have addressed the role of M cells, 
located in the midmyocardium. The prolonged action potential of M cells leads to dispersion of repolarization 
which is increased significantly in patients with the long QT syndrome. We have focused on the long QT 
3 syndrome which is marked by an increase in the late Inh and have investigated two previously proposed 
mechanisms, using numerical simulations of inhomogeneous tissue cables and sheets. 

The first mechanism proposes that early after depolarizations in the M cell domain leads to premature 
stimuli in surrounding tissue. However, our numerical results show that under normal electrotonic coupling, the 
late sodium current is not sufficiently strong to induce an excitation in neighboring points. Typical currents 
that are responsible for the generation of early after depolarizations are fairly small. 3 For example, in the 
case studied here, the late I^a current is on the order of 2-4 /iAcm~ 2 while in long QT 1 and long QT 2, the 
respective reduction of Ik s and Ikt is on the order of 1 /iAcm~ 2 . Furthermore, the L type calcium current 
involved in the generation of early after depolarizations is roughly of the same magnitude. This should be 
compared to the strength of the injection current needed to induce a counter propagating wave in a cable with 
fi = 0. Numerically, we find that the necessary injection current is approximately an order of magnitude larger 
(<~ 30 /iAcm~ 2 ). It is therefore perhaps not a surprise that the inclusion of a leak current in I^a is not sufficient 
to create a wave propagating in the unidirectional direction. We have verified that a reduction of Ikt was also 
unable to initiate unidirectional propagation. 

We also investigated the effect of a discontinuity in conductivity between the M cell domain and the epi- 
cardium, as reported in a study. 8 We found that a large increase in resistance between the two domains could 
result in an extra stimulus in the epicardium that was provoked by an early after depolarization in the M cell 
domain. However, for the model used in our study, this increase was at least four times higher that the one 
reported in the literature. This indicates that the discontinuity in conductivity between the M cell domain and 
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epicardium is likely insufficient to generate a premature stimulus arising from the M cell domain. 

Our simulation study focused on early after depolarizations originating in M cells. We did not consider other 
possible sources of early after depolarizations. For example, cells in the right ventricular epicardium exhibit 
a prominent transient outward current /to- 35 Early after depolarizations originating from the right ventricle 
epicardial cells might play an role in phase 2 reentry, 36 which is thought to be associated with the Brugada 
syndrome. 37 A recent numerical study, using a modified version of the LRd model that included I to together 
with a large reduction in the time constant for the L type calcium channel, investigated this scenario. 38 These 
modifications were found, under very specific conditions, to elicit a secondary wave in a cable divided into two 
sub-domains with different /to densities. However, this scenario is specific to right ventricular epicardial cells 
and is unlikely to play a role in reentry phenomena associated with M cell domains in the left ventricle. 

The second mechanism we investigated relies on the slow repolarization of the M cell domain which functions 
as a wave block for subsequent stimuli. When the M cell domain is inhomogeneous, or alternatively, if the pacing 
breaks the translational symmetry along the endocardial wall, a properly timed extra stimulus can be blocked in 
one area of the myocardium while allowed to propagate in another. This, then, could create reentrant excitation. 

We first studied the propagation of stimuli, delivered prematurely with varying coupling intervals, for dif- 
ferent basic cycle length and long QT syndrome strength. We found a rich set of possible dynamics of the 
wave following the premature stimulus. Particularly noteworthy was the dramatic reduction of the wave speed 
occurring within the M cell domain. Prerequisite for this behavior is a sufficiently long basic cycle length. 
Preliminary investigations of the dynamics indicate that during this slow propagation the usual sodium channel 
activation of neighboring cells fails. Instead, the calcium channel becomes responsible for the depolarization 
and the spread of the impulse. This scenario has also been observed experimentally in growth cultures of rat 
myocytes. 39 

The ID simulations revealed two separate regions in the basic cycle length vs. /i space in which conduc- 
tion block could occur. The first region, characterized by relatively small basic cycle length, exhibits a small 
conduction block window. For coupling intervals within this window, the wave fails to propagate abruptly. 
The second region, characterized by a slow wave propagation and a subsequent wave block, has a much larger 
conduction block window. Furthermore, along the boundary of this region, the dynamics is not straightforward. 
A window of coupling intervals that exhibit conduction block is followed by a window of coupling intervals that 
show slow propagation within the M cell domain. Again, the calcium handling seems to be a likely candidate 
for an explanation of this behavior. This is supported by recent experimental work, able to visualize voltage 
and calcium simultaneously, that showed that calcium dynamics plays a crucial role in the generation of early 
after depolarizations in particular and action potential prolongation in general. 40 Thus, a model that faithfully 
incorporates calcium handling is essential. This might require a description of local spatial events including 
calcium sparks, 41 currently absent from the LRd model. 

In mapping out the phase diagram of Fig. [5] we have focused on a point in the parameter space within this 
second region (basic cycle length of 1000 ms and fi=0.08). The main finding of our 2D simulations is that the 
second mechanism can lead to reentrant excitations. Indeed, for large enough w, or large enough r, we find 
reentry that can last 1 s or more. The observed pseudo-ECG, with its polymorphic undulating morphology, is 
suggestive of torsade de pointes. Furthermore, for a large bubble size, the reentry becomes pinned to the region 
of heterogeneity, thus leading to long lasting reentrant excitations. 

We should note that the observed mechanism for reentrant excitations is not specific to long QT 3. We have 
also investigated both the first and second mechanism in a model for long QT 2. In this model, we reduced the 
maximal conductance of Ikt which allowed us to vary the strength of the long QT syndrome. As in our long QT 
3 model, we found that early after depolarizations did occur in the M cell domain but were not able to excite 
neighboring tissue. On the other hand, when an M cell domain was included in the transmural sheet, the long 
QT 2 model was able to partially block a premature stimulus and sustained reentrant excitations were found. 
Also, we point out that we have replicated the pacing protocol employed in Ref. 11. There, the endocardium 
was paced at a long basic cycle length, followed by a premature stimulus on the epicardium. As the pacing 
protocol already breaks the translational symmetry along the endocardium, we included a continuous strip in 
our numerical geometry. We found that this pacing protocol could lead to reentry, provided that the strip had a 
sufficient width. Besides the geometrical construction of this paper or the pacing protocol of, 11 there are other 
ways to break the translational symmetry. For example, a strip of M cells with domains of different long QT 
syndrome strength can have, for a given coupling interval, domains that exhibit wave block while other domains 
support normal propagation. This, again, would result in partial wave block and subsequent reentry. 
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The mechanism that was shown here to create reentrant excitation relies on an extra stimulus, delivered at a 
relatively small coupling interval. This extra stimulus can possibly be generated by an early after depolarization 
in a Purkinje fiber that successfully excites the endocardium. A possible mechanism that would lead to the 
initiation of a wave originating from a Purkinje fiber early after depolarization is based on the experimentally 
observed reduced coupling between the Purkinje fiber and the endocardium. 42 Under this condition, numerical 
simulations have revealed that the inclusion of a large resistive barrier between the Purkinje fibers and the 
ventricular cells can lead to early after depolarizations in the Purkinje fibers that are able to generate extra 
stimuli in the ventricular domain. 31 

Our numerical simulations were carried out using a simplified transmural geometry of the ventricular wall 
that consisted of a centrally located M cell domain, sandwiched between two electrophysiologically equal domains 
representing the endo- and epicardium. Incorporating more realistic geometries and endo- and epicardial cells 
with different electrophysiological properties, however, will most likely not change the qualitative conclusions 
of this study. Furthermore, the action potentials of the epicardial and endocardial cells recorded in the wedge 
preparation of the canine left ventricle, were very similar. 11 

The numerical simulations employed 2D transmural sheets while the actual heart is of course a much more 
complicated 3D object. Nevertheless, our numerical work, replicating the experiments in wedges, 11 demonstrates 
that the creation of an intramural reentry event is plausible under certain conditions. In 3D, an ensuing 
intramural filament can either close on itself, leading to a scroll ring that remains intramural and is on the 
scale of the ventricles, 43,44 or can reorient and attach itself to the endo- and epicardium. Possible subsequent 
instabilities could destabilize this filament, resulting in multiple filaments and fibrillation. Only full-scale 3D 
simulations, preferably on realistic geometries containing detailed anatomical information, can address these 
events. We are currently investigating this possibility, which includes the development of novel algorithms that 
can accurately simulate electrical wave propagation in anatomically realistic hearts. 45 
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